
#include "LinearWaveDGKernel.h"
#include "ReconstructionUserObject.h"

registerMooseObject("O4App", LinearWaveDGKernel);

template <>
InputParameters
validParams<LinearWaveDGKernel>()
{
  InputParameters params = validParams<DGKernel>();
  return params;
}

LinearWaveDGKernel::LinearWaveDGKernel(const InputParameters & parameters)
  : DGKernel(parameters),
	VFVInterface(*getCheckedPointerParam<VFVProblem*>("_vfv_problem"), _assembly)
{

}

void LinearWaveDGKernel::computeReconstructionResidual(Moose::DGResidualType type)
{
	mooseInfo("LinearWaveDGKernel::computeReconstructionResidual");
}

void LinearWaveDGKernel::computeReconstructionJacobian(Moose::DGJacobianType type)
{
	mooseInfo("LinearWaveDGKernel::computeReconstructionJacobian");
}

Real LinearWaveDGKernel::computeQpResidual(Moose::DGResidualType type)
{
  Real r = 0;

  RealGradient _velocity = RealGradient(1,0,0);
  Real vdotn = _velocity * _normals[_qp];

  Real ul = _u_rs[_qp];
  Real ur = _u_rs_neighbor[_qp];
  Real flux = ul * vdotn;
  switch (type)
  {
    case Moose::Element:
		r = flux * _test[_i][_qp];
      break;

    case Moose::Neighbor:
        r = -flux * _test_neighbor[_i][_qp];
      break;
  }

  return r;
}

Real LinearWaveDGKernel::computeQpJacobian(Moose::DGJacobianType type)
{
   Real r = 0;
   RealGradient _velocity = RealGradient(1,0,0);

   Real vdotn = _velocity * _normals[_qp];

  switch (type)
  {
    case Moose::ElementElement:
      if (vdotn >= 0)
        r += vdotn * _phi[_j][_qp] * _test[_i][_qp];
      break;

    case Moose::ElementNeighbor:
      if (vdotn < 0)
        r += vdotn * _phi_neighbor[_j][_qp] * _test[_i][_qp];
      break;

    case Moose::NeighborElement:
      if (vdotn >= 0)
        r -= vdotn * _phi[_j][_qp] * _test_neighbor[_i][_qp];
      break;

    case Moose::NeighborNeighbor:
      if (vdotn < 0)
        r -= vdotn * _phi_neighbor[_j][_qp] * _test_neighbor[_i][_qp];
      break;
  }

  return 0;
  //return r;
}
